      PROGRAM FLO83

!     ******************************************************************
!     *                                                                *
!     *   PROGRAM TO SOLVE THE TWO DIMENSIONAL EULER EQUATIONS         *
!     *                                                                *
!     ******************************************************************

!     ******************************************************************
!     *                                                                *
!     *   ANALYSIS OF TRANSONIC FLOW PAST AN AIRFOIL ON AN C-MESH      *
!     *                                                                *
!     ******************************************************************

!     ******************************************************************
!     *                                                                *
!     *   COPYRIGHT (C) ANTONY JAMESON 1988                            *
!     *                                                                *
!     ******************************************************************

!     ******************************************************************
!     *                                                                *
!     *   CELL-CENTERED MULTIGRID SCHEME                               *
!     *                                                                *
!     ******************************************************************

!     ******************************************************************
!     *                                                                *
!     *   REFERENCES :                                                 *
!     *   (1) A. JAMESON,"SOLUTION OF THE EULER EQUATIONS BY A         *
!     *       MULTIGRID METHOD",APPLIED MATHEMATICS AND COMPUTATION,   *
!     *       VOL. 13,1983,PP. 327-356                                 *
!     *   (2) A. JAMESON,"TRANSONIC FLOW CALCULATIONS FOR AIRCRAFT",   *
!     *       IN NUMERICAL METHODS IN FLUID DYNAMICS,EDITED BY         *
!     *       F. BREZZI,LECTURE NOTES IN MATHEMATICS,VOL. 1127,        *
!     *       SPRINGER VERLAG,1985,PP.156-242                          *
!     *   (3) A. JAMESON,"MULTIGRID ALGORITHMS FOR COMPRESSIBLE FLOW   *
!     *       CALCULATIONS",PROC. SECOND EUROPEAN CONFERENCE ON        *
!     *       MULTIGRID METHODS,COLOGNE,1985,EDITED BY U. TROTTENBERG  *
!     *       AND W. HACKBUSCH,LECTURE NOTES IN MATHEMATICS,VOL. 1228, *
!     *       SPRINGER VERLAG,1986,PP. 166-201                         *
!     *   (4) A. JAMESON,"COMPUTATIONAL TRANSONICS",                   *
!     *       COMM. PURE APPL. MATH.,VOL. 16,1988,PP. 507-549          *
!     *   (5) A. JAMESON,"COMPUTATIONAL AERODYNAMICS FOR AIRCRAFT      *
!     *       DESIGN",SCIENCE,VOL. 245,1989,PP. 361-371                *
!     *                                                                *
!     ******************************************************************

!     PROGRAMMED BY ANTONY JAMESON,SUMMER 1987
!     SUCCESSOR TO FLO52,PROGRAMMED BY ANTONY JAMESON ,1980-1985

!     THE SOLUTION IS CALCULATED ON A SEQUENCE OF MESHES FOR M=1,MMESH
!     USING A MULTIGRID CYCLE FOR EACH MESH

!     ******************************************************************
!     *                                                                *
!     *   THE PROGRAM HAS FIVE MAIN MODULES                            *
!     *                                                                *
!     *               INPUT AND CONTROL MODULE                         *
!     *               MESH GENERATOR                                   *
!     *               FLOW SOLVER                                      *
!     *               MULTIGRID MODULE                                 *
!     *               OUTPUT MODULE                                    *
!     *                                                                *
!     *   THE FLOW SOLVER AND MULTIGRID MODULE ARE CONFIDENTIAL        *
!     *                                                                *
!     ******************************************************************

!     ******************************************************************
!     *                                                                *
!     *   COMMENTS IN SUBROUTINE INPUT DESCRIBE THE INPUT PARAMETERS   *
!     *   AND SPECIFY RECOMMENDED VALUES FOR GENERAL USE               *
!     *                                                                *
!     ******************************************************************
      USE CARBUNCLE
      IMPLICIT NONE

!     ******************************************************************
!     *                                                                *
!     *   FLOW VARIABLES                                               *
!     *                                                                *
!     ******************************************************************

!     W(I,J,1)    = DENSITY
!     W(I,J,2)    = MOMENTUM IN X DIRECTION
!     W(I,J,3)    = MOMENTUM IN Y DIRECTION
!     W(I,J,4)    = TOTAL ENERGY

!     ******************************************************************
!     *                                                                *
!     *   MULTIGRID CONTROL PARAMETERS                                 *
!     *                                                                *
!     ******************************************************************

!     NMESH         THE NUMBER OF MESHES IN THE MULTIGRID CYCLE
!     LCYC.EQ.1     V CYCLE
!     LCYC.EQ.2     W CYCLE
!     KXPAND.EQ.1   DESCEND TO A COARSER MESH
!     KXPAND.EQ.-1  ASCEND TO A FINER MESH
!     MODE.EQ.0     THE FINEST MESH IN THE CYCLE
!     MODE.EQ.1     ANY COARSER MESH IN THE CYCLE
!     KODE.EQ.1     PERFORM A NORMAL UPDATE OF THE SOLUTION
!     KODE.EQ.2     FIRST ENTRY TO A COARSER MESH DURING THE CYCLE
!     KODE.EQ.-1    CALCULATE RESIDUALS WITHOUT UPDATING THE SOLUTION
!     FCOLL         RELAXATION FACTOR FOR RESIDUAL COLLECTION
!     FADD          SMOOTHING FACTOR FOR MULTIGRID CORRECTIONS
!     FBC.EQ.1.     UPDATE THE FAR FIELD ON THE COARSE MESHES
!     FBC.EQ.0.     FREEZE THE FAR FIELD ON THE COARSE MESHES

!     ******************************************************************
!     *                                                                *
!     *   DIMENSION CONTROL PARAMETERS                                 *
!     *                                                                *
!     ******************************************************************

!      INTEGER, PARAMETER :: IDM=385, JDM=65                  ! MOVED TO MODULE CARBUNCLE
!      INTEGER, PARAMETER :: IDN=IDM+1, JDN=JDM+1             ! MOVED TO MODULE CARBUNCLE
!      INTEGER, PARAMETER :: IDP=IDM+2, JDP=JDM+2             ! MOVED TO MODULE CARBUNCLE
!      INTEGER, PARAMETER :: IDQ=(IDM-1)/2+2, JDQ=(JDM-1)/2+2 ! MOVED TO MODULE CARBUNCLE
!      INTEGER, PARAMETER :: IMX=385                          ! MOVED TO MODULE CARBUNCLE

!     ******************************************************************
!     *                                                                *
!     *   MULTIGRID STORAGE CONTROL PARAMETER                          *
!     *                                                                *
!     ******************************************************************

!     192*32 MESH :     IDX = 2500
!     384*64 MESH :     IDX = 9096

!      INTEGER, PARAMETER :: IDX=9096, IDX2=2*IDX, IDX4=4*IDX    ! MOVED TO MODULE CARBUNCLE

!     ******************************************************************
!      COMMON/DAT/ GAMMA,RM,RHO0,P0,EI0,H0,C0,U0,V0,CA,SA,ALPHA          ! MOVED TO MODULE CARBUNCLE
!      COMMON/LIM/ NX,NY,IL,JL,IE,JE,ITL,ITU                             ! MOVED TO MODULE CARBUNCLE
!      COMMON/PRM/ CFL,VT,HM,VIS0,VIS2,VIS4,ADIS,BC,                     ! MOVED TO MODULE CARBUNCLE
!     .            SMOOPI,SMOOPJ,KSMOOP,MSTAGE
!      COMMON/OUT/ CP(IDM),SCAL,CHORD,XM,YM,CL,CD,CM                     ! MOVED TO MODULE CARBUNCLE
!      COMMON/CNV/ RTRMS,RTMAX,IRT,JRT,HRMS,HMAX,IH,JH,NSUP              ! MOVED TO MODULE CARBUNCLE
!      COMMON/CYC/ MCYC,NCYC,NTIM,NOUT,NPRNT,IPRNT,LPRNT,NMESH,K1,KX,KW  ! MOVED TO MODULE CARBUNCLE
!      COMMON/MGR/ KODE,MODE,FCOLL,FADD,FBC                              ! MOVED TO MODULE CARBUNCLE
!      COMMON/MGQ/ CFLF,CFLC,HMF,HMC,LCYC                                ! MOVED TO MODULE CARBUNCLE
!      COMMON/MON/ TOT,RTRMS0,NRES,IMON                                  ! MOVED TO MODULE CARBUNCLE
!      COMMON/STT/ ISTART                                                ! MOVED TO MODULE CARBUNCLE
!      COMMON/SEQ/ FCYC(4),FPRNT(4),FOUT(4),FTIM(4),                     ! MOVED TO MODULE CARBUNCLE
!     .            GPRNT(4),HPRNT(4),HMESH(4),MMESH
!      COMMON/TIT/ TITLE(20)                                             ! MOVED TO MODULE CARBUNCLE

      INTEGER :: KCYC(7),IWRIT,M
      REAL :: TIM

!     ******************************************************************

      IWRIT     = 6
      TIM       = CUTIME(0)
      WRITE (IWRIT,700) TIM
      WRITE (IWRIT,600)
      WRITE (IWRIT,2)
    2 FORMAT(14H0PROGRAM FLO83,30X,&
             35HANTONY JAMESON,PRINCETON UNIVERSITY/&
             35H0AIRFOIL ANALYSIS IN TRANSONIC FLOW,&
             40H BY SOLUTION OF UNSTEADY EULER EQUATIONS)

!     ******************************************************************
!     *                                                                *
!     *   SET UP THE AERODYNAMIC AND GEOMETRIC DATA                    *
!     *                                                                *
!     ******************************************************************

      CALL INPUT
      K1        = 1
      KX        = 1
      KW        = 1
      MODE      = -1

!     INCREASE THE DIMENSIONS TO THOSE OF THE FINEST MESH
!     IN THE SEQUENCE

      IF (MMESH.EQ.1) THEN
         CALL GMESH       ! GENERATE THE MESH FOR THE FINAL CALCULATION
      ELSE
        DO 10 M=2,MMESH
           NX = NX + NX
           NY = NY + NY
   10   CONTINUE
        CALL GMESH        ! GENERATE THE MESH FOR THE FINAL CALCULATION
      END IF

!     ******************************************************************
!     *                                                                *
!     *   GENERATE THE MESH FOR THE FINAL CALCULATION                  *
!     *                                                                *
!     ******************************************************************

!   11 CALL GMESH    ! COMMENTED OUT - INCLUDED ABOVE

!     EXPAND THE MESH TO PRODUCE THE FIRST MESH IN THE SEQUENCE

      IF (MMESH.EQ.1) THEN
         M = 1
      ELSE
        DO 20 M=2,MMESH
           CALL MOVCO
   20   CONTINUE
        M = 1
      END IF

!     SET UP THE METRIC AND INITIAL DATA

   31 CALL SETUP
  start from here***********************************************************
C     ******************************************************************
C     *                                                                *
C     *   CALCULATE THE FLOW                                           *
C     *                                                                *
C     ******************************************************************
C
      TIM       = CUTIME(0)
      WRITE (IWRIT,700) TIM
      MCYC      = FCYC(M)
      NTIM      = FTIM(M)
      NOUT      = FOUT(M)
      NPRNT     = FPRNT(M)
      IPRNT     = GPRNT(M)
      LPRNT     = HPRNT(M)
      LPRNT     = MAX(1,LPRNT)
      NMESH     = HMESH(M)
      DO 40 N=1,NMESH
      KCYC(N)   = 0
   40 CONTINUE
      L         = (MCYC  -1)/500  +1
      NOUT      = MAX(NOUT,L)
      NRES      = 0
      NCYC      = 0
      TOT       = 0.
      GCL       = FCL
      NSTART    = 0
      IF (RM.GT.1..AND.M.EQ.1) NSTART = 100
      N         = NMESH
      IMON      = 0
C
C     PRINT DATA TO MONITOR CONVERGENCE
C
      CALL MONITR
      IMON      = 1
C
C     ******************************************************************
C     *                                                                *
C     *   START OF THE MULTIGRID CYCLE                                 *
C     *                                                                *
C     ******************************************************************
C
  101 NCYC      = NCYC  +1
      ISTART    = 0
      IF (NCYC.LE.NSTART) ISTART = 1
      KODE      = 1
      MODE      = 0
      CFL       = CFLF
      HM        = HMF
      WORK      = 1.
C     NSUP      = -1
C     IF (NCYC.EQ.1.OR.MOD(NCYC,NPRNT).EQ.0) NSUP = 0
      NSUP      = 0
  111 KXPAND    = 1
      KCYC(N)   = KCYC(N)  +1
  121 TOT       = TOT  +WORK
C
C     CALCULATE THE PERMISSIBLE TIME STEP
C
      IF (MOD(NCYC-1,NTIM).EQ.0) CALL STEP
C
C     UPDATE THE FLOW
C
      FCL       = 0.
      IF (MOD(NCYC-1,5).EQ.0) FCL = GCL
      CALL EULER
      KODE      = 1
      IF (N.LT.NMESH.OR.MOD(NCYC-1,NOUT).NE.0) GO TO 131
C
C     PRINT DATA TO MONITOR CONVERGENCE
C
      CALL MONITR
  131 IF (NCYC.GE.MCYC) GO TO 201
      IF (MOD(NCYC,NPRNT).EQ.0.AND.N.EQ.NMESH) CALL OUTPUT
      IF (KXPAND.EQ.-1.AND.N.LT.NMESH) GO TO 141
      IF (N.EQ.1) GO TO 101
      KODE      = -1
C
C     CALCULATE THE RESIDUALS
C
      CALL EULER
C
C     TRANSFER TO THE NEXT COARSER MESH
C
      CALL MOVCO
      KODE      = 2
      MODE      = 1
      N         = N  -1
      CFL       = CFLC
      HM        = HMC
      WORK      = .25*WORK
      IF (N.GT.1) GO TO 111
      KXPAND    = -1
      GO TO 121
C
C     TRANSFER TO THE NEXT FINER MESH
C
  141 CALL MOVFIN
      N         = N  +1
      WORK      = 4.*WORK
      IF (N.EQ.NMESH) GO TO 101
      IF (KCYC(N).LT.LCYC) GO TO 111
      KCYC(N)   = 0
      GO TO 141
C
C     ******************************************************************
C     *                                                                *
C     *   END OF THE MULTIGRID CYCLE                                   *
C     *                                                                *
C     ******************************************************************
C
C     PRINT DATA TO MONITOR CONVERGENCE
C
  201 IMON      = 2
      CALL MONITR
C
C     ******************************************************************
C     *                                                                *
C     *   POSTPROCESS THE RESULT                                       *
C     *                                                                *
C     ******************************************************************
C
      CALL OUTPUT
      TIM       = CUTIME(0)
      WRITE (IWRIT,700) TIM
C
C     ******************************************************************
C     *                                                                *
C     *   TRANSFER TO A FINER MESH IF M.LT.MMESH                       *
C     *                                                                *
C     ******************************************************************
C
      IF (M.GE.MMESH) STOP
      M         = M  +1
      CALL MOVFIN
      GO TO 31
C
C     ******************************************************************
C
  600 FORMAT(1H1)
  610 FORMAT(1X,6F12.4)
  630 FORMAT(1H0,17A4)
  640 FORMAT(1X,6I12)
  700 FORMAT(1H0,5HTIME ,F8.2,7HSECONDS)
      END
